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Abstract 

Binary supermassive black holes form naturally in galaxy mergers, but their long-term evolution is uncertain. 
In spherical galaxies, A/-body simulations show that binary evolution stalls at separations much too large for 
significant emission of gravitational waves (the "final parsec problem")- Here, we follow the long-term evolu- 
tion of a massive binary in more realistic, triaxial and rotating galaxy models. We find that the binary does not 
stall. The binary hardening rates that we observe are sufficient to allow complete coalescence of binary SBHs 
in 10 Gyr or less, even in the absence of collisional loss-cone refilling or gas-dynamical torques, thus providing 
a potential solution to the final parsec problem. 
Subject headings: 



1. INTRODUCTION 

When two galaxies containing supermassive black holes 
(SBHs) merge, a binary SBH forms at the center of the new 
galaxy. The two SBHs can eventually coalesce, but only 
after stellar- or gas-dynamical processes bring them close 
enough together ( <, 10~ 2 pc) that gravitational radiation is 
emitted. There is strong circumstantial evidence that rapid 
coalescence is the norm. For insta nce, no binar y SBH has 
ever b een unambiguously observed (Merritt & Milosavljevic 
120051) . Furthermore, in a galaxy containing an uncoalesced 
binary, mergers would eventually bring a third SBH into 
the nucleus, precipitating a gravitational slingshot interac- 
tion th at would eject one or mor e of the SB Hs from the nu- 
cleus dMikkola & Valtonenll990tlVolonteri. Haardt & Madaul 
120031) . This could produce off-center SBHs, and could 
also weaken the tight correlations that are observed between 
SBH mass and galaxy properties (Ferrarese & Merritt 2000; 
Graham et al.l200lHMarconi & Huni2003l) . 

Unless the binary mass ratio is extreme, dynamical friction 
rapidly brings the smaller SBH into a distance ~ G/Lt/o 2 from 
the larger SBH, where /j = M\Mz/{M\ +M2) is the binary re- 
duced mass and a is the ID velocity dispersion of the stars. 
At this separation - of order 1 pc - the two SBHs begin to 
act like a "hard" binary, ejecting passing stars with velocities 
large enou gh to remove them from the nucleus. jV-body sim- 
ulations dMakino & Funatol 120041 ISzell Merritt & M ikkola 
120051 iBerczik. M erritt & Spurzem 120051) show that contin- 
ued hardening of the binary takes place at a rate that depends 
strongly on the number N of "star" particles used in the simu- 
lation. As N increases, the hardening rate falls, as expected if 
the binary's loss cone is re populate d by star-star grav itational 
encounters (Yu 2002; Milosa vlievic & Merrittll2003l) . When 
extrapolated to the much larger N of real galaxies, these re- 
sults suggest that binary evolution would generally stall (the 
"final parsec problem"). 

To date, A/-body simulations of the long-term evolution of 
binary SBHs have only been carried out using spherical or 
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nearly spherical galax y models. But it has been suggested 
(Merritt & Poon 2004) that binary hardening might be much 
more efficient in non-axisymmetric galaxies due to the quali- 
tatively different character of the stellar orbits. Here, we test 
that suggestion by carrying out the first A^-body simulations of 
massive binaries in strongly non-axisymmetric galaxy mod- 
els. We find that the hardening rate is independent of N for 
particle numbers up to at least 0.4 x 10 6 . To the extent that 
our galaxy models are similar to real merger remnants, these 
results imply that binary SBHs can efficiently harden through 
purely stellar-dynamical interactions in many galaxies, thus 
providing a plausible solution to the final parsec problem. 

2. METHOD 

Our iV-body models were generated from the phase-space 
distribution function 

f(E,L z ) = const x (e~^ E - l) e" pn °^ (1) 

( iLaeoute & Longarettilll996l) . Here E = v 2 /2 + 4> is the en- 
ergy per unit mass of a star, z) is the gravitational poten- 
tial in the meridional plane, and L z is the angular momentum 
per unit mass in the direction of the symmetry (z) axis; the 
potential is set to zero at the radius where the density falls to 
zero. The quantity in parentheses on the right h and side of 
equation ([0 is the energv-dependent lKing l lfl966) distribution 
function. The additional, angular-momentum-dependent fac- 
tor has the effect of flattening the models and simultaneously 
giving them a net rotation about the z axis. The degree of flat- 
tening can be specified via the dimensionless rotation param- 
eter ©0 = y/9 j (4nGpo)Qo, with po the central mass density 
of the galaxy. The parameter P determines the central concen- 
tration of the model; its value was chosen such that the spher- 
ical isotropic model generated from equ ation iffll ha d a dimen- 
sionless central concentration Wq = 6 l lKingll966l) . Here and 
below we adopt standard A^-body units, i.e. the gravitational 
constant and total mass of the galaxy are one, and the galaxy's 
energy is —1/4. 

A pair of massive particles representing the two SBHs were 
introduced into the models at time t = 0. The two particles 
were given equal masses, M\ = Mi = M./2, and were placed 
on coplanar, circular orbits at distances ±0.3 from the galaxy 
center in the equatorial plane. In most of the simulations de- 
scribed below, M. = 0.0 4. This is rather larger th an the typ- 
ical ratio, ~ 1 x 10~ 3 ( Merri tt"& FerraresdlZOOll) . observed 
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between SBH mass and galaxy mass; such a large mass for 
the SBH particles was chosen in order to minimize the rate 
of relaxation-driven loss cone refilling, which occurs mor e 
rapidly for smaller M. dBerczik. Merritt & Spurzeml 120051) . 
and to come as close as possible to the "empty loss cone" 
regime that characterizes real (axisymmetric) galaxies. In or- 
der to estimate the dependence of the binary decay rate on M., 
we carried out a limited set of additional simulations with dif- 
ferent values of M., as described below. 

Integrations of the particle equations of motion were 
carried out using a high-accuracy, d irect-summation A^-body 
code (Berczik. Merritt & Spurzem 2005) on two paral- 
lel sup ercomp uters incorpor at ing sp ecial-pupose GRAPE 
jFukushige. Makino & KawaH 120051) accelerator boards: 
gravitySimulator 5 and GRAC E. 6 Integration parameters 
were s imilar to those adopted in B erczik. Merritt & Spurzeml 
(2005|and we refer the reader to that paper for details about 
the performance of the code. Integrations were carried out 
for various N in the range 0.025 <iV<0.4x 10 6 and for 
various values of the galaxy rotation parameter coo in the 
range < COo < 1.8. In addition, each model was integrated 
with two choices for the orientation of the binary's angular 
momentum, either parallel to that of the galaxy ("prograde") 
or counter to it ("retrograde"). 

3. RESULTS 

After the two SBH particles come close enough together 
to form a bound pair, the parameters of their relative Ke- 
plerian orbit can be computed. Figure la shows the evo- 
lution of I /a, the binary inverse semi-major axis, in a set 
of simulations with COo = and various N. These spheri- 
cal, non-rotatin g models are very similar to th e models con- 
sidered in Berczik. Merritt & Spurzem (2005), and the bi- 
nary evolution found here exhibits the same strong N depen- 
dence that was observed in that study: the hardening rate, 
s(t) = (d/dt)(\/a), is approximately constant with time and 
decreases rough ly as N ~ 1 . This behavior has been described 
quantitatively (Milosavlievic & Merritt 2003) on the basis of 
loss-cone theory: stars ejected by the binary are replaced in a 
time that scales as the two-body relaxation time, and the latter 
increases roughly as N in a galaxy of fixed mass and size. 

When the rotation parameter COo is increased to ~ 0.6, the 
initially axisymmetric models become unstable to the forma- 
tion of a bar, yielding a slowly-tumbling, triaxial spheroid. 
Figure 2 illustrates the instability via snapshots of the COo = 
1.8 model integration. This model is moderately flattened 
initially, with mean short-to-long axis ratio of ~ 0.46, and 
strongly rotating, with roughly 40% of the total kinetic energy 
in the form of streaming motion. Movies based on the simula- 
tion reveal that the two SBH particles initially come together 
by falling inward along the bar before forming a bound pair. 

The long-term behavior of the binary (Figure lb) is strik- 
ingly different in this rotating model than in the spherical 
model: not only is the hardening rate high, but more signif- 
icantly, it shows no systematic dependence on particle num- 
ber. In fact, the two simulations of the COo = 1.8 model with 
largest N (200fc and 400£) exhibit almost identical evolution 
of the binary. 

To determine the dependence of the binary's evolution rate 
on the properties of the galaxy model, a suite of A^-body in- 

5 http://www.cs.rit.edu/ grapecluster/clusterlnfo/grapeClusterlnfo.shtml 
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FIG. 1. — Evolution of the binary inverse semi-major axis, \/a, in N- 
body simulations with various N. (a) Spherical, nonrotating galaxy model 
((Bo = 0). (b) Flattened, rotating galaxy model (COo = 1.8). At f w 10, this 
model forms a triaxial bar (cf. Figure 2). 
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FIG. 2. — Snapshots at four times (/ = 0,5, 10,30) of the particle positions, 
projected onto the (x.y) plane, in an /V-body integration with COo — 1-8 and 
N = 200k. The SBH particles are indicated in green. The evolution of the 
binary semi-major axis in this integration is shown in Figure lb. 



tegrations were carried out for various (coo,A0 and for M. = 
0.04. Axis ratios of the galaxy were compu ted using its mo- 
ment o f inertia tensor, as described in Dub inski & Carlbergl 
( 1199 11) . The results are summarized in Figure 3. After a 
strong bar forms at f w 10 in the unstable models, it evolves 
gradually toward rounder shapes, but the system maintains 
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time N 
FIG. 3. — (a) Evolution of the intermediate-to-long axis ratio b/a of the galaxy in W-body integrations with N = 200k and various values of the rotation 
parameter ct>o; Mi = M2 = 0.02. For (Do > 0.6, the galaxy is unstable to non-axisymmetric deformations and forms a slowly-tumbling triaxial spheroid. The 
length c of the short axis is nearly constant with time and is determined by the initial flattening of the model, i.e. by (Do. The thin black line shows the result of an 
integration that omitted the two SBH particles, (b) Hardening rate of the binary measured at t = 150. The value of (Do is indicated via the key on the lower left; 
black symbols are "prograde" integrations, i.e. the binary revolves in the same sense as the galaxy, and red symbols are "retrograde" integrations. The dashed 
lines shows the N~ 1 dependence that is characteristic of a collisionally-resupplied loss cone. 



a significant triaxiality until the end of the simulation. This 
sl ow evolution appears sim ilar to that in the simulations 
of Theis & Spurzem ( 1999), where two-body relaxation was 
identified as the driving mechanism. The presence of a mas- 
sive binary in our s imulatio ns might also tend to d estroy 
the triaxiality ( Athanassoula, Lambert & Dehnen 2005), al- 
though an integration excluding the binary showed a similar 
degree of evolution; see Figure 3a. 

The steep, ~ A^ 1 dependence of the binary hardening rate 
in the spherical model changes to an essentially constant hard- 
ening rate for COo > 1.2 (Figure 3b). Even the ©0 = 0.6 
model - which, Figure 3a suggests, is close to marginal sta- 
bility - yielded a hardening rate that was substantially larger 
than in the sph erical models, consistent with suggestions 
jMerritt & Poonl 120041) that even slight departures from ax- 
isymmetry could significantly influence the binary's evolu- 
tion. The hardening rate was found not to depend on the initial 
sense (prograde vs. retrograde) of the binary orbit (Figure 3b). 

Since M,/M ga i = 0.04 is considerably larger than the ra- 
tio ~ 1 x 10~ 3 observed in re al galaxies JMerritt & Ferraresel 
120011 iMarconi & Huntl 12003). we carried out an additional 
set of simulations in order to evaluate the M, -dependence of 
the binary hardening rate. These integrations used ©0 = 1.8, 
N = 0.1 x 10 6 , and 0.01 < M. < 0.08. We found that the 
hardening rate increased with decreasing M„. At t = 100, 
the hardening rates were s = (20.0, 13.0,8.2,3.4) for M. = 
(0.01,0.02,0.04,0.08). These results should be interpreted 
with caution since we did not vary N and therefore can not 
state with certainty whether the s values are ^-independent. 
The linear size of the binary's loss cone is proportional to M. 
making it easier for two-body scattering to affect the harden- 
ing rate as M. is decreased. But the M. = 0.04 integrations 
are clearly not in the collisionally-repopulated regime (Figure 
lb) and so the differences that we observe between the hard- 
ening rates for M. = 0.04 and 0.08 are likely to be robust. 
Based on these results, it is reasonable to conclude that binary 



evolution would be substantially more rapid than implied by 
Figure 3b in the case M. w 10~ M ga i. 

4. DYNAMICAL INTERPRETATION 

A hardening rate that is independent of N implies a colli- 
sionless, i.e. relaxation-independent, mode of loss-cone re- 
filling. Just such a mode is expected in triaxial galaxies: 
the lack of an axis of symmetry implies that stellar orbits 
need not conserve any component of the angular momentum, 
hence they can pass arbitrarily close to t he center after a finite 
time and interact with a central object (tNorman & SiUdlT983l 
Gerhard & Binney 1985). 

A full derivation of the expected rate of supply of stars to 
the binary in these models would be very difficult, but we can 
do an approximate calculation. The rate per unit of orbital 
energy at which centrophilic orbits supply mass (i.e. stars) to 
a region of radius r t at the center of a galaxy is 

M(E)dE = r,A(E)M c (E)dE (2) 

where A(E)d is the rate at which a single star on a centrophilic 
(e.g. box or chaotic) orbit of energy E experiences near-center 
passages with pericenter distances < d, and M c (E)dE the 
mass in stars on centrophilic o rbits with energies from E to 
E + dE jMerritt & Poonl2004l) . Setting r, = Ka, with K w 1, 
gives the mass flux into t he binary's sphere of influence; the 
implied hardening rate is (Berczik. Merritt & Spurzem 2005) 

dt \a J aM, J 

Here, (C) « 1.25 is the average value of the dimensionless 
energy change during a single star-binary encounter, C = 
[M./2m*](AE/£). 
Our Af-body models have density p ~ r~ 2 beyond the core 
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radius r c w 0.25. In a p « /• galaxy, 

A (E) » ^ e -(^-£/.)/- 2 , Mc {E)=f c (E) x ^^ e (^-^)/2- 2 
»* 9 G 

(4) 

with r/, = GM./o 2 , £), = 4>(r/,), and frfg) the fraction of th e 
orbits at energy E that are centrophilic ( Merr itt & Poonl2004l) . 
The implied binary hardening rate is 



4^6 {C)Kf c 



on 



e -(E-EH)/2J dE ~ 2 .5f£. (5) 



Here f c is an energy-weighted, mean fraction of centrophilic 
orbits, and the lower integration limit was set to E„; the latter 
can only be approximate since the true density of our galaxy 
models departs from r~ 2 at r < r c w r/,. Substituting a w 0.47 
and r/, s» 0.18 from the galaxy models gives s as 40/ c . By 
comparison, the hardening rates in the A^-body models reach 
a peak value at t w 20 of s w 16, consistent with the de- 
rived expression if f c w 1 /2. The gradual drop observed in 
the hardening rate at later times, s(t) w 16. — 5.21n(f/20), 
20 < f < 250, suggests that the number of centrophilic stars 
is becoming smaller, due to depletion by the binary and to the 
gradual change in the galaxy's shape. 

Taken at face value, equation Q implies s M, 2 ; however 
for small M., r/, <C r c and the assumption that p ~ r~ 2 , r > r/, 
breaks down. In any case, the observed dependence of s on 
M, is slightly weaker, s ~ M~ l (§3). 

5. IMPLICATIONS 

The time s cale for gravi tational wave emission by a binary 
black hole is (Peters 196*1) 



\6 A F{e) o 8 M 2 



a h 



(6) 



Here M. = Mi +M2, /J = MyM^jM 1 is the reduced mass of 
the binary, a is the ID central velocity dispersion of stars 
in the nucleus, and a/, = G/j/4o 2 is the semi-major axis of 
the binary when it first becomes "hard," i.e. tightly bound; 
the factor F(e) depends on the binary's orbital eccentricity 
andF(0) = 1. In order that gravitational-wave-driven coales- 
cence take place in less than 10 10 yr, an equal-mass, circular- 
orbit binary wit h M. = 10 8 Mm must first rea ch a separa- 
tion a <, 0.05a/, (iMerritt & Milosavlie'vic1l2005l) . This is just 
achieved in our simulations: at, w 1.1 x 10~ 2 , and the final 
value of a in the bar-unstable models is ~ 6 x 10~ 4 . This 



is a conservative interpretation since (1) for reasonable scal- 
ings of our galaxy model to real galaxies (e.g. total mass 
= 10 n M o , half- mass radius = 10 3 pc), an elapsed time of 
250 in A^-body units corresponds to < 1 Gyr; (2) the binary 
is continuing to harden at the final time-step in our simula- 
tions (Figure lb); (3) our experiments with different M. found 
s ~ M~ 1 , implying substantially more rapid hardening in the 
case M,/M ga i w 10~ 3 ; (4) the binary had nonzero eccentricity 
in our simulations. In addition, gas is a significant component 
of disk galaxies and, in many mergers, the final hardening 
of the binary would be acc elerated by gas-dynamical torques 
( Escalaetal 12001 IPottTet alJ2005l) . 

Our simulations of binary evolution are substantially more 
realistic than existing ones based on spherical or nearly- 
spherical galaxy models. Even more realistic simulations, 
which follow both the early and late stages of a merger 
between two galaxies, are probably beyond the capabili- 
ties of current algorithms and hardware due to the need to 
accurately treat both large (~ 10 kpc) and small (~ 0.01 
pc) spatial scales. However, our galaxy models (slowly- 
tumbling triaxial spheroids ) are similar to those produced 
in full merger simulations JBournaud. Jog & Combesl 120051 
Naab, Khochfar & Burkert 2006), suggesting that our results 
for the long-term evolution of the binary are probably fairly 
generic in spite of the rather artificial initial conditions. 

Uncertainties about the resolution of the "final parsec prob- 
lem" have been a major impediment to predicting the fre- 
quency of SBH mergers in galactic nuclei, and hence to com- 
puting event rates for proposed gravitational wave interfer- 
ometers like LISA. 8 If binary coalescence rates are assumed 
to be similar to galaxy merger rates, gravitational wave events 
integrated o ver the o bservabl e universe could b e as frequent 
as 10 2 yr" 1 <Haehnelll99"llSesana et ail20 04). Our results, 
combined with the indirect evidence that binary SBH coales- 
cence is efficient, suggest that such high event rates should be 
taken seriously. 
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